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We report an ensemble nuclear magnetic resonance (NMR) implementation of a quantum lattice 

O ' 



o 

(N 



gas algorithm for the diffusion equation. The algorithm employs an array of quantum information 
processors sharing classical information, a novel architecture referred to as a type-11 quantum com- 
puter. This concrete implementation provides a test example from which to probe the strengths 
and limitations of this new computation paradigm. The NMR experiment consists of encoding a 
mass density onto an array of 16 two-qubit quantum information processors and then following 
the computation through 7 time steps of the algorithm. The results show good agreement with 
the analytic solution for diffusive dynamics. We also describe numerical simulations of the NMR 
implementation. The simulations aid in determining sources of experimental errors, and they help 
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■ define the limits of the implementation. 
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I. INTRODUCTION 
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The advent of fast quantum algorithms |l| has spawned a broad search for new algorithms that utilize the novel 
features of quantum information. Among the new proposals are quantum lattice gas (QLG) algorithms, which, in 
analogy to their classical counterparts, make use of arrays of interacting sites to perform useful calculations. In 
the quantum case, however, the sites behave quantum mechanically, while the site-to-site interactions can be either 

. New algorithms have been devised to solve selected computational problems 
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2] or quantum mechanical 3] 
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classical 

such as the diffusion equation 



In the case where the quantum mechanical sites 



and the Schrodingcr equation 
(or nodes) communicate with each other classically, the required architecture for QLG algorithms has been termed a 
type-II quantum computer^ |. 

A tvp e-II device is essentially a classically parallel computer, with the exception that the computing elements 
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follow the rules of quantum mechanics. The advantage gained from the classical network is completely analogous 
to the improvement gained in a classical, massively-parallel architecture. However, the use of quantum mechanical 
nodes introduces several notable differences. Classical lattice gas algorithms become unstable (and unusable) when 
the relevant transport coefficient is reduced or when nonlinearities are increased. In the quantum case, the transport 
coefficient and degree of nonlinearity can be varied at will using the appropriate quantum operations at each site. In 
addition, the quantum algorithms typically require a smaller number of qubits per site than do the classical algorithms. 
Finally, the family of QLG algorithms can handle increasingly complex calculations as the number of (qu)bits per 
site is increased. For example, when two qubits are present in a site, the QLG algorithms can solve the relatively 



simple diffusion equation in one, two, or three dimensions 
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4 and the more difficult nonlinear Burgers equation in one 



dimension 25]. With four qubits per site, the QLG algorithms can solve coupled nonlinear field equations governing 
(he vrli icily and magnetic fields of one-dimensional magnetohydrodynamic turbulence 2f|. With six qubits per site, 
the QLG algorithms can model the nonlinear Navier-Stokes equations in two dimensions governing a viscous fluid ^]. 
A more complete description of type-II quantum computers and their scaling properties has been given by Jeffrey 



re compl' 
Yepez |25j|2g. 



Here, we present a methodology for implementing a quantum lattice gas algorithm on a nuclear magnetic resonance 
(NMR) 9] type-II architecture. In this implementation, we encode a discrete mass density onto distinct spatial 
locations of a liquid-state sample. We use magnetic field gradients to discriminate between locations in the sample, 
thus creating an array of addressable ensemble NMR quantum information processors. In addition, we use radio 
frequency (RF) pulses and methods learned from previous work 3 3 to execute the necessary operations in each 
quantum processor. The result is a concrete implementation examining the necessary control for realizing a quantum 
lattice-gas algorithm using NMR techniques. 



II. LATTICE GAS ALGORITHMS 

The lattice gas method is a tool of computational physics used to model hydrodynamical flows that are too large for 
a standard low-level molecular dynamics treatment and that contain discontinuous interfacial boundaries that prevent 



a high-level partial differential equations description 



The basic idea underlying the lattice gas method 



is to statistically represent a macroscopic scale time-dependent field quantities by "averaging" over repeated instances 
of a system of artificial microscopic particles scattering and propagating throughout a lattice of interconnected sites. 
A particular instance of the system has many particles distributed over the lattice sites. Multiple particles may coexist 



at each site at a given time, and each particle carries a unit mass and a unit momentum of energy. Particles interact 
on site via an artificial collision rule which exactly conserves the total mass, momentum, and energy at that site. 
The movement of particles along the lattice is prescribed by a streaming operation that shifts particles to nearest 
neighboring sites, thus endowing the particles with the property of momentum. The lattice gas algorithm encapsulates 
the microscopic scale kinematics of the particles scattering on site and moving along the lattice. The mean-free path 
length between collisions is about one lattice cell size and the mean-free time between collision elapses after a single 
update. This is computationally simple in comparison to molecular dynamics where many thousands of updates are 
required to capture such particle interactions. 

The mesoscopic evolution is obtained by taking the ensemble average over many instances of microscopic realizations. 



At the mesoscopic scale, the average presence of each particle type is defined by an occupation probability. In addition, 
the microscopic collision and streaming rules translate into the language of kinetic theory 16, 
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The behavior 



of the system is described by a transport equation for the occupation probabilities, and this equation is a discrete 



Boltzmann equation called the lattice Boltzmann equation li 
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The lattice Boltzmann equation further translates into a macroscopic, continuous, effective field theory by letting 
the cell size approach zero (the limit of infinite lattice resolution called the continuum limit). At the macroscopic 
scale, partial differential equations describe the evolution of the field, admitting solutions such as propagating sound 
wave modes and diffusive modes. The passage of the Boltzmann equation to the effective field theory begins by 
expanding the occupation probabilities, which have a well-defined statistical functional form, in terms of the continuous 
macroscopic variables, such as the mass density p (and the velocity or energy field if they are defined in the model). 
This expansion usually is carried out perturbatively in a small parameter such as the Knudsen number (ratio of 
mean-free path to the largest characteristic length scale) or the Mach number (ratio of the sound speed to the largest 
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characteristic flow speed) in a fashion analogous to the Chapman-Enskog expansion of kinetic theory 
Conversely, and self-consistently, the macroscopic field quantities can also be expressed as a function of the mesoscopic 
occupation probabilities-for example, the mass density at some point is a sum over the occupation probabilities in 
that vicinity. 

Quantum lattice gas algorithms are generalizations of the classical lattice gas algorithms described above where 
quantum bits are used to encode the occupation probabilities and where the principle of quantum mechanical super- 
position is added to the artificial microscopic world. In this quantum case, the mesoscopic occupation probabilities are 
mapped onto the wave functions of quantum mechanical sites. In the case where the quantum lattice gas describes a 



hydrodynamic system when the time evolution of the flow field is required, we must periodically measure these occu- 
pation probabilities and the quantum lattice gas algorithm becomes suited to a type-II implementation. Such type-II 
algorithms have been shown to solve dynamical equations such as the diffusion equation the Burgers equation 
1 2 J , and magnetohydrodynamic Burgers turbulence ^(j . As a first exploration of a type-II architecture using NMR, 
we implemented a QLG model of diffusive dynamics in one dimension. 

III. SOLVING THE 1-D DIFFUSION EQUATION 

The quantum lattice gas algorithm that solves the 1-D diffusion equation derives from a classical lattice gas of 

. n 

particles moving up and down a 1-D lattice|J|. The motion of the particles occurs in discrete steps (streaming phase), 
and the particles have a probability of changing directions (collision). When the collisions are such that the particles 
reverse directions half of the time, then the continuum effective field theory that emerges obeys diffusive dynamics. 
In this case, the motion of an individual particle is a random walk, and an arbitrary initial distribution of particles 
will diffuse isotropically as a function of time. 

The lattice gas described above is summarized by the Boltzmann equation 

f h2 (z ± Az,t + At) = f h2 (z,t) + n h2 (z,t), (1) 

where the left-hand side denotes the occupation of the lattice as a function of the previous lattice configuration and 
where the collision term is 



fil,2=±i[/l(l-/2)-/ 2 (l-/l)] (2) 



The variables f 1 = f 1 (z,t) and f 2 = f 2 (z,t) are the occupation probabilities for finding upward- and downward- 
moving particles, respectively, at the site location z and time t. The time step is denoted by At, while the lattice 
spacing is given by Az. The collision term changes the direction of some particles, and it is responsible for the diffusive 
behavior. 

The interesting macroscopic quantity of the lattice gas is the mass density field, p, defined as the sum of upward- 
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and downward-moving particles 

p(z,t) = h(z,t) + f 2 (z,t). (3) 

The ambiguity in assigning the mass density between the two occupation probabilities is resolved by a constraint for 
local equilibrium demanding that the mass density be initially distributed equally 

f?(n, 0) = fp(n, 0) = \ P {nAz, 0). (4) 

After a single time step, the occupation probabilities fx and /2 evolve according to l[T]l. resulting in a new mass density 

p(z, t + At) = i [p(z + Az, t) + p(z - Az, t)} . (5) 

The first finite-difference in time of the mass density field is then written as 

p(z, t + At)- p(z, t) = i [p(z + Az, t) - 2p{z, t) + p(z - Az, t)} . (6) 

In the limit where the lattice cell size and the time step approach zero (Az — > and At — * 0), the mass density field 
becomes continuous and differentiable. The second-order Taylor expansion of equation (jSJ) about z and t can thus be 
written in the differential form 



dp{z,t) _ (Az 2 \ d 2 p(z,t) 



dt V 2At / dz 2 

where it is now evident that p evolves according to the diffusion equation with a constant transport coefficient ^Ei- 
Finally, in this implementation we consider an initial mass density p(z, t = 0) whose evolution obeys the periodic 
boundary condition p(z,t) = p(z + L,t), where L is the length of the lattice. As a result, the initial mass density 
diffuses until the total mass is evenly dispersed throughout the lattice. 

The corresponding quantum lattice gas algorithm description begins by encoding the occupation probabilities, and 
thus the mass density, in the states of a lattice of quantum objects. The streaming and collision operations are 
then a combination of classical and quantum operations, including measurements. The aim of the algorithm is to 
take an initial mass density field and to evolve its underlying occupation probabilities according to the Boltzmann 
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equation Q). A schematic of the entire quantum algorithm is shown in Fig. 1. A single time step of the algorithm is 
decomposed into four sequential operations: 

1 . encoding of the mass density 

2. application of the collision operator C at all sites 

3. measurement of the occupation numbers 

4. streaming to neighboring sites. 

These operations are repeated until the mass density field has evolved for the desired number of time steps. In the 
first time step, the encoding operation specifies the initial mass density profile, while in all the subsequent steps the 
encoding writes the results of the previous streaming operation. The final time step ends with the readout of the 
desired result, so operation 4 is not performed. 

Each occupation probability is represented as the quantum mechanical expectation value of finding a two-level 
system, or qubit, in its excited state |1). As a result, the state of the qubit encoding the value f a (z, t) is 

\fa(z,t)) = Vfa(z,t)\l) + y/l- f a (z,t)\0). (8) 

It follows that a single value of the mass density is recorded in two qubits, one for each occupation number. The 
combined two-qubit wave function for a single node becomes 

M*,*)> = v / 7i^|ii) + v / /i(i-/ 2 )|io>+ (9) 
V(i-/i)/ 2 |oi) + x/(i-/i)(i-/ 2 )|oo). 

The kets 1 00) , 1 01) , |10), and 1 11) span the joint Hilbert space of the two qubits, and this is the largest dimension space 
over which quantum superpositions are allowed. As with the classical algorithm, the constraint for local equilibrium 
i jlj l forces the initial occupation probabilities at a node to be half of the corresponding mass density value. 

The occupation numbers encoded in the two-qubit wave function \ip(z,t)) can be recovered by measuring the 
expectation value of the number operator h a , as given in 



fa(z,t) = {1>(z,t)\n a \il>(z,t)), 



(10) 



where hi = h ® 1, hi = 1 §5 h, where 1 is the 2x2 identity matrix, and where the action of the single-qubit number 
operator h returns 1 if qubit is in its excited state and for the ground state. 

The encoded occupation probabilities evolve as specified by the Boltzmann equation by the combined action of the 
collision operator, the measurement, and streaming. The collision operator contributes by taking the local average 
of the two occupation probabilities. This averaging (not to be confused with statistical coarse-grain averaging, time 
averaging, or ensemble averaging) is done by choosing the the collision operator C to be the "square-root of swap" 
gate, written as 



C = 



1 
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(11) 



in the standard basis. The same collision is applied simultaneously at every site, resulting in 



\ij'(z,t))=C\i>(z,t)). 



(12) 



Using l(TU|) . the intermediate occupation probabilities of the wave function \ip'(z,t)) are 



(13) 



as required for a = 1,2. The third operation physically measures these intermediate occupation probabilities f' a (z,t) 
at all the sites. If the algorithm is performed on individual quantum systems, then the values are obtained by 
averaging over many strong quantum measurements of identical instances of each step. However, when the algorithm 
is performed using a sufficiently large ensemble of quantum systems, as in the case of NMR, then a single weak 
measurement of the entire ensemble can provide sufficient precision to obtain f' a {z, t). A single time step is completed 
with the streaming of the occupation probabilities to the nearest neighbors, according to the rule 



f 1 (z-Az,t + At)=fi(z,t) 
f 2 (z + Az,t + At) = &(z,t) 



(14) 
(15) 
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The information on the two qubits is shifted to the neighboring sites in opposite directions. The streaming operation 
is a classical step causing global data shifting, and it is carried out in a classical computer interfaced to the quantum 
processors. Together, the last three operations result in 

/t,a(* ±Az,t + At) = i [f ± (z, t) + f 2 (z, t)] 

(16) 

which is the exact dynamics described by the Boltzmann equation (|T)). 

IV. NMR IMPLEMENTATION 
A. Spin System and Control 

The goal of the NMR implementation is to experimentally explore the steps outlined by the diffusion QLG algorithm. 
For this two-qubit problem, we chose a room-temperature solution of isotopically-labeled chloroform ( 13 CHCi3), where 
the hydrogen nucleus and the labeled carbon nucleus served as qubits 1 and 2, respectively |27j. The chloroform sample 
was divided into 16 classically-connected sites of two qubits each, creating an accessible Hilbert space larger than 
would be available with 32 non-interacting qubits. 

The internal Hamiltonian of this system in a strong and homogeneous magnetic field Bq is 

Haemal = ~ {IhBq) (lcB Q ) <j\ + (17) 

where the first two terms represent the Zeeman couplings of the spins with B and the last term is the scalar coupling 
between the two spins. The operators of the form are Pauli spin operators for the spin a and the Cartesian 
direction k. The choice of chloroform is particularly convenient because the different gyromagnetic ratios, jh and 
7c, generate widely spaced resonant frequencies. As a result, a RF pulse applied on resonance with one of the spins 
does not rotate, to a very good approximation, the other spin. In the 7 T magnet utilized for the implementation, the 
hydrogen and carbon frequencies were about 300 MHz and 75 MHz, respectively. The widely spaced frequencies allow 
us to write the two RF control Hamiltonians as acting on the two spins independently. The externally-controlled RF 
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Hamiltonians are written as 



H a RF {t) = -\ [w a x (t)a a x + w a y (t)a a y ] . (18) 



The RF Hamiltonians generate arbitrary single-spin rotations with high fidelity when the total nutation frequencies 



sai sci, 

33 
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are much stronger than J, the scalar coupling constant. The scalar coupling Hamiltonian and the single-spin rotations 
permit the implementation of a universal set of gates, and they are the building blocks for constructing more involved 
gates such as the collision operator C 

The lattice of quantum information processors is realized by superimposing a linear magnetic field gradient on the 
main field Bo, adding a position dependent term to the Hamiltonian having the form 

tt ( \ 1 ( dB z \ , I ( dB z X 2 

H gradlent {z) = -- I lH-g^z \a z - - [Jc-q^z I <V (20) 

The variable z denotes the spatial location along the direction of the main field, while the constant specifies the 
strength of the gradient. The usefulness of this Hamiltonian can be appreciated by noticing that the offset frequencies 
A£Ih,c — 1h,c (^§r) z °f the spins vary with position when the gradient field is applied. Spins at distinct locations 
can thus be addressed with RF fields oscillating at the corresponding frequencies. In this way, the magnetic field 
gradient allows the entire spin ensemble to be sliced into a lattice of smaller, individually addressable sub-ensembles. 

Using the coupling, RF, and gradient Hamiltonians described above, together with the appropriate measurement 
and processing tools, we can now describe in detail how the four steps of the diffusion QLG algorithm translate to 
experimental tasks. The lattice initialization step (1) uses the magnetic field gradients to establish sub-ensembles 
of varying resonant frequency addressable with the RF Hamiltonians. The collision step (2 ) m akes use of both the 



step (2 ) m 

n cliabj 



The readout (3) is 



RF and the internal coupling Hamiltonians to generate the desired unitary operatior 
accomplished by measuring the spins in the presence of a magnetic field gradient. And finally, the streaming operation 
(4) is performed as a processing step in a classical computer in conjunction with the next initialization step. 
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B. Lattice Initialization 

The initialization of the lattice begins by transforming the equilibrium state of the ensemble into a starting state 
amenable for quantum computation. At thermal equilibrium, the density matrix is 



&the 



exp 



B- internal 
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(21) 



where e has a value on the order of 10 -5 . The equilibrium state is highly mixed and the two spins have unequal 
magnetizations. To perform quantum computations, it is convenient to transform the equilibrium state into a pseudo- 
pure state 2^, 2{]], a mixed state whose deviation part transforms identically to the corresponding pure state and, when 
measured, returns expectation values proportional to those that would be obtained by measuring the underlying pure 
state. Two transformations create the starting pseudo-pure state |00) from the thermal state. First, the magnetizations 
of the two spins are equalized, 



^thermal 



Equalize 



1 e 



2 2 2 



Vequal = 7T7 + 77 ( 1 + — ) Wl + a 



1C 



(22) 



followed by a pseudo-pure state creation sequence that results in 



Pseudo-pure 1 , v/3 

&equal &pp ^ 



v» r- I - 1 \ a z + a z + a z G z 

■2- 1\ 2 V ic. ' 



(23) 



The equalization and pseudo-pure state creation sequences are described in detail in reference [30j. For clarity, we 
define the constant in front of the brackets to be e', allowing us to write the pseudopure state a pp in terms of the 
desired spinor 1 00) as 



1 



e' l + e'|00)(00|. 



(24) 



Expressed in this manner, it is now seen that a unitary transformation applied to a pp acts trivially on the term 
proportional to the identity, but it evolves the term |00)(00| as it would a pure state. 

Individually addressing the sites of the lattice, as depicted in Fig. 1, is accomplished by selectively addressing 
p „ f *e samP , TU P ^ , ^ ,o sM ec t ,» ta _ „ sg ,„ g 

(MRI) 31], and it works by applying the gradient Hamiltonian in the presence of suitably shaped RF pulses. First, 
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consider the Hamiltonian for a one-spin system subjected to a linear magnetic field gradient in the ^-direction and to 
a time-dependent RF pulse applied in the y-direction. In this case, the Hamiltonian is 



H 



RF,G 



Cm) 



1 ( dB z 

az 



-w y (t)a v 



(25) 



where the a z term is the linearly-varying static field and the a y term is the time-dependent RF. The Hamiltonian 
HnF.G( z T~t) does not commute with itself at all times, so a closed- form and exact solution cannot be easily given 
without specifying the function w y (t). A valuable approach, however, is to consider the approximate evolution 
generated by Hnp j Q(z,t) during infinitesimal periods of the RF pulse. To first order, the evolution during the initial 
period At becomes 



URF,c{z,t = At) « exp 



1 / dB z A 

l 2 r— A ' • 



exp 



.Wy(At)At 



(26) 



By defining the term in the parenthesis as Ak z = ^y^^-At, the evolution of an initial density matrix a z through a 
single period becomes 



Urf^zU^pq ps exp 



Ak z 



a x exp 



Ak z 



w y (At)At + a z 



(27) 



where small angle approximations have been made. The first term is a spatial helix of the x and y magnetizations 
having a wavenumber Ak z . The second term is the first order approximation to the magnetization remaining in the 
state a z . Another period of evolution will affect the a z term as described, creating a new magnetization helix with 
wavenumber Ak z . In addition, the initial helix will have its wavenumber increased by an amount Ak z . The final 
result over many periods is the formation a shaped magnetization profile having many components 



N 



nAk z z 

l z cr 2 



er^exp 



nAk z z 
-i — ^—o- 2 



w y (nAt)At 



(28) 



Each term in summation can be interpreted as a cylindrical Fourier component of the x-y magnetization weighted by 

the RF nutation rate w v {nAt). The RF waveform specifies the magnitude of each spatial Fourier component, and the 

l I 

resulting spatial profile is the Fourier transform of the RF waveform 32J. An equivalent description is to say that, for 
weak RF pulses, the excited magnetization of the spins at a given resonance frequency is, to first order, proportional 
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to the Fourier component of the RF waveform at that frequency. As a result, control of the appropriate RF Fourier 
component essentially translates to selective addressing of spatial frequencies, which in turn allows the excitation of 
particular spatial locations. 

The Fourier transform approximation allows encoding of arbitrary shapes on the various spatial locations of one 
uncoupled nuclear species. For QIP, however, coupled spins are required to implement two-spin operations. In 
particular, the chloroform carbons and protons are coupled together via the scalar coupling. Given that the required 
RF waveforms should be weak, the coupling interferes with the desired evolution. The effect of the coupling present 
while encoding on spin 1 is removed by applying a strong RF decoupling sequence on the second spin [34l |. The 
decoupling modulates the cr z operator in the interaction Hamiltonian, making its average over a cycle period equal 
to zero. As a result, the second spin feels an identity operation during the decoupling. Fig. 2 shows the complete 
RF and gradient pulse sequence. As can be seen from the diagram, the first encoding on qubit 1 was subsequently 
swapped to qubit 2, followed by a re-encoding of qubit 1. We chose this method because the smaller gyromagnetic 
ratio of 13 C causes a narrower frequency dispersion in the presence of the gradients, making the carbon decoupling 
simpler. 

As described above, the encoding process writes the desired shapes in the spatial dependence of each spin's x- 
magnetization. The occupation numbers, however, are proportional to the z-magnetization, as can be seen when the 
number operator in the equation 



f a (n,m) = (ip(n,m)\h a \ip(n 7 m)), 



(29) 



is replaced with n, 



a — 



i(l + cr°) resulting in 



f a (n,m) = - [1 + (^(n,m)K|V(n,m))] . 



(30) 



where second term in the brackets represents the z-magnetization. The encoding process is followed by a tt/2 pulse 
that rotates the excited x-magnetization to the z direction. 
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C. Collision and Swap Gates 

After initialization, the next step is to apply the collision operator. For the QLG algorithm solution to the diffusion 
equation, the collision operator C is the square-root of swap gate. Expressed in terms of the Pauli operators, it is 



C = exp 



• f 1 2 12 1 2^ 



(31) 



where an irrelevant global phase has been ignored. Written in this form, the operation C can be decomposed into a 



sequence of implementable RF pulses and scalar coupling evolutions 
the exponent commute with each other, resulting in 



ll| by noticing that the product operators in 



C — exp 
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(32) 



Expanding the first and last exponentials as scalar couplings sandwiched by the appropriate single-spin rotations 
results in 



C = 



exp 
exp 
exp 
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(33) 
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The exponents of terms proportional to cr^er! represent internal Hamiltonian evolutions lasting for a time t°°* = 1/(4 J). 
The exponents of terms with single-spin operators are implemented by ir/2 rotations. They were generated by RF 
pulses whose nutation rate was about 50 times greater than J. All of the pulses and delays were applied without a 
magnetic field gradient in order to transform all of the sites identically. 

As shown in Fig. 2, swap gates were utilized both in the lattice initialization and in the measurement of the carbon 
magnetization. The pulse sequence for the swap gates was almost identical to the sequence for C. The only difference 
was that the internal evolution delay was set to t s z ™ ap = 1/(2 J). 
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D. Measurement 

The occupation numbers resulting from the collision were obtained by measuring the z-magnetizations and using 
equation i|30|) . Since only the <r" and a® operators are directly observable, a "read out" w/2 pulse transformed the 
z-magnetization into x-magnetization. The proton magnetization was measured directly after the collision, while the 
carbon magnetization was first swapped to the protons before observation. Measurements of both the 13 C and 1 H 
magnetizations were carried out separately, and in both cases via the more sensitive proton channel. The measurements 
were made in the presence of a weak linear magnetic field gradient, causing signals from different sites to resonate with 
distinguishable frequencies. The observed proton signal was digitized and Fourier transformed to record an image of 
the spatial variation of the spin magnetization. The observed spectrum was then processed to correct the baseline and 
to obtain the resulting magnetization at each site. Because each site is composed of a slice of the sample with spins 
resonating in a band of frequencies, the occupation number for each site was obtained by averaging over all spins in 
the corresponding band. 

E. Streaming 

The final step involves classically streaming the results of the measurements according to equations Ijl4|l and l|15(l . 
The streaming operation is applied in conjunction with the next lattice initialization step by adding a linearly varying 
phase to the Fourier transform of the desired shape. The added phase causes a shift in the frequency of the pulse 
determined by the slope of the phase. When the frequency-shifted pulse is applied in the presence of the magnetic field 
gradient, the shift results in spatial translation of the encoded shape. The streaming operation is thus implemented 
as a signal processing step in the lattice initialization procedure. 

V. RESULTS AND DISCUSSION 

The results of the experiment are shown in Fig. 3, together with plots of the analytical solution and of numerical 
simulations of the NMR experiment. In total, 7 steps of the algorithm were completed using a parallel array of 16 
two-qubit ensemble NMR quantum processors. The observed deviations between the data points and the analytical 
plots can be attributed to imperfections in the various parts of the NMR implementation. 

To explore the source and relative size of these errors, we simulated perfect experiments, each time adding controlled 
errors in four sections of the implementation: 
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• Fourier transform approximation in the initialization 

• Decoupling during the initialization 

• Encoding swap gate and 7r/2 pulse errors 

• Collision gate errors 

The Fourier transform approximation executes a correct writing of the desired magnetization to first order in the 
overall flip angle. To explore errors introduced by the approximation, we simulated NMR experiments using nutation 
angles ranging from 7r/2 to 7r/20. In this range, angles smaller than 7r/4 resulted in accurate encodings of the desired 
Gaussian shapes through the ten steps of the implementation. The errors in the three remaining sections were 
simulated by using RF pulses with the actual time and nutation rate that were used on the spectrometer. By using 
a finite power, errors from imperfect averaging of the scalar coupling could be observed. Errors in the collision gate 
caused the least impact to the mass density, followed by errors originating from the imperfect decoupling sequence. 
The largest deviations originated from realistic simulations of the swap gate and the 7r/2 pulses in the encoding. It 
is important to note that the simulated gate fidelities for the swap and collision gates, although imperfect, are still 
about 0.995. This suggests that the observed deviations are caused by the coherent buildup of errors through a few 
iterations, and not just by the individual errors from a single gate. The complete simulation, using realistic RF pulses 
and a shaped pulse nutation angle of 7r/4, is plotted in Fig. 3. The calculated mass densities closely match the 
experimental results, suggesting that the observed errors are accurately modeled. 

Other potential sources of errors include the finite signal to noise, the state fidelity of the starting pseudo pure 
state, and gradient switching time. In addition, spin relaxation, random self-diffusion of the liquid molecules, and RF 
inhomogeneity can all cause attenuations in the strength of the signal. In our experiments, these last three mechanisms 
manifested themselves indirectly through reduced signal to noise. However, given that this attenuation was present in 
all of the experiments, any direct results were mostly normalized away in the data processing. Although none of the 
above errors contributed significantly to our implementation, they arc likely to become important as more complicated 
algorithms are executed on larger lattices. 

In particular, molecular diffusion over the time of an operation places a lower bound on the physical size of the 
volume element corresponding to each site in the computation. In the 1-D case discussed here, the root-mean-squared 
displacement (Az = \j2Dt) for chloroform (D = 2.35 x 10~ 5 cto 2 /s) is about 10.8^m over the 25ms needed for 
encoding and the collision operator. Since the actual volume element were about 625^m across, this resulted in a 
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negligible mixing of the information in adjacent sites. However, it is clear that for this approach to type-II quantum 
computer to remain viable for large matries and more complex collision operators the physical size of the sample must 
grow with the size of the problem. 



VI. CONCLUSION 

Ensemble NMR techniques have been used to study the experimental details involved in quantum information 
processing. The astronomical number of individual quantum systems (~ 10 18 ) present in typical liquid-state spin 
_ bleS ^ ^ to. ^ „ f _ g spto u ^ ,e ™ ble mtu „ ta 

been successfully utilized to create the necessary pseudopure states |2F 



operations over the ensemble 




22( and to systematically generate nonunitary 



In this experiment, we again exploit the ensemble nature, but this time as a means 
of realizing a parallel array of quantum information processors. The novel architecture is then used to run a quantum 
lattice gas algorithm that solves the 1-D diffusion equation. 

The closeness of the data to the analytical results is encouraging, and it demonstrates the possibility of combining 
the advantages of quantum computation at each node with massively parallel classical computation throughout the 
lattice. Currently, commercial MM machines routinely take images with 256 x 256 x 256 volume elements. As a 
result, the large size of the NMR ensemble provides, in principle, sufficient room to explore much larger lattices. 
However, in moving to implementations with more computational power, several challenges remain. The limited 
control employed here is sufficient for a few time steps of the algorithm, but refinements are necessary to increase 
the number of achievable iterations. In addition, although complicated operations have been done in up to 7 NMR 
qubitsj^, 3^, 37 1, the problem of efficiently initializing a large lattice of few-qubit processors still remains. Our results 
provide a first advance in this direction, and they provide confirmation that NMR techniques can be used to test 
these new ideas. 
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FIG. 1: The circuit diagram shows the quantum lattice gas algorithm for solving the 1-D diffusion equation. The algorithm 
employs N two-qubit sites to encode the discretized mass density. Each site codes for a single value of the mass density 
using the quantum state of the two qubits. The encoded information is subjected to a series of local transformations that 
evolve the system. The collision operator C is the only potentially entangling operation in the algorithm, and it creates 
quantum coherences limited to each two-qubit system. The streaming is executed by classical communication, and it moves 
the occupation numbers up and down the lattice as denoted by the arrows. The sectioned cylinder depicts the position of the 
adjacent sites in the NMR sample. Each site is physically realized as an addressable slice of isotopically-labeled Chloroform 
solution. 
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FIG. 2: The NMR implementation consists of four main sections, each corresponding to the prescribed QLG algorithm step. 
The top two lines in the diagram correspond to RF pulses applied to the proton and carbon qubits, respectively. The third 
line shows the application of magnetic field gradients. In the encoding section, the initial carbon magnetization is recorded on 
the protons before being transferred to the carbons. The starting magnetization is specified by using a RF pulse shaped as the 
Fourier transform of the desired magnetization. The shaped pulses are applied in the presence of gradients so that each site can 
be addressed. A carbon decoupling sequence prevents the scalar coupling from interfering with the low power shaped pulses. 
The 7r/2 at the end of the encoding move the information form the x-axis to the z-axis, as required by the QLG algorithm. The 
collision operator follows the encoding, and it is implemented without gradients to ensure that all of the sites in the sample 
feel the same transformation. The results are observed in two experiments, each time using the more sensitive proton channel. 
A swap gate is added when measuring the carbon magnetization. Finally, the streaming operation is applied by shifting the 
frequencies of the carbon and proton shapes in opposite directions. 
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FIG. 3: The experimental mass densities are plotted in the figure, together with plots of the analytical solution and the 
numerical simulation of the NMR experiment. The normalized, dimensionless mass densities are plotted as they were encoded 
on the lattice. Seven steps of the algorithm were implemented on 16 two-qubit sites. The simulations were performed using 
the actual RF nutation rates and times of the experimental setup. The calculations closely match the data, suggesting that 
the deviation between the analytical results and the data can be attributed imperfections in the methodology. As a result, the 
simulations promise to be useful in exploring the errors from alternate methods. 



